Spread of correlations in long-range interacting systems 
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Understanding the dynamics of many-body systems is crucial for understanding, e.g., thermaliza- 
tion or transmission of information. Nevertheless, little is known in the case of quantum systems 
with long-range interactions. Here, we analyze the long-range Ising model in a transverse field, 
where interactions decay as a power-law with distance oc r~^, a > 0. Using complementary nu- 
merical and analytical techniques, we identify three dynamical regimes: short-range-like with an 
emerging light cone for a > 2; weakly long-range for 1 < a < 2 without a clear light cone but 
with a finite propagation speed of excitations; and fully non-local for a < 1 with instantaneous 
transmission of correlations. This last regime breaks generalized Lieb-Robinson bounds. Numerical 
calculation of the entanglement spectrum demonstrates that the usual picture of propagating quasi- 
particles remains valid for long-range interactions. This allows an intuitive interpretation in terms 
of qualitative changes to the spin-wave dispersion, leading to diverging quasi-particle velocities in 
the long-range regime. Our results may be tested in state-of-the-art trapped-ion experiments. 



Recent years have seen a surge of interest in the dy- 
namics of quantum systems. One reason is that clean 
quantum dynamics can now be observed in finely con- 
trolled physical setups, such as ultracold atoms [T]-[3]. 
Other experiments in this direction can be expected in 
the near future in polar molecules [IHS] and Rydberg 
atoms [3 [8], or with trapped ions [9UT4]. where first steps 
towards studying dynamics have already been done [M]. 
It is currently not clear, however, how the long-range in- 
teractions in these cases influence the dynamics. 

From the theoretical side, one seeks to understand the 
long-time behavior of out-of-equilibrium many-body sys- 
tems, especially driven by the question how thermaliza- 
tion can emerge from unitary quantum dynamics. The 
hope is that, after long times, a small region would no- 
tice the presence of the rest of the system as if it was in 
contact with a thermal bath and would essentially equi- 
librate to a thermal state. We already know that not all 
systems equilibrate to thermal states [15], which imme- 
diately calls for a more general classification of possible 
equilibrium states. For example, it is now widely ac- 
cepted that integrable models equilibrate locally towards 
Generalized Gibbs states [16 . 

The typical scenario for such studies are global 
quenches., where an initial state is evolved with an Hamil- 
tonian it is not an eigenstate of. Unfortunately, apart 
from rare cases where models are exactly solvable, for 
global quenches we do not yet have good tools to address 
the long-time regime with classical computers O [TTl [18] . 
All tools that have successfully been applied to static 
properties of quantum many-body systems face specific 
difficulties when dealing with out-of-equilibrium dynam- 
ics. Results obtained with Monte Carlo in imaginary time 
need to be continued analytically to real time, which is 



by no means easy and sometimes even impossible due 
to dynamical quantum phase transitions [19, 20 ; sim- 
ulations based on Tensor Network (TN) states such as 
Matrix Product States (MPS), or the Multi-Scale En- 
tanglement renormalization Ansatz (MERA), typically 
require a small amount of entanglement to be accurate, 
while during time evolution correlations build up linearly 
in time as a consequence of Lieb-Robinson bounds [21] , 
and soon overcome the limit where the numerical simu- 
lations are feasible [22 . 

A simpler class of out-of-equilibrium dynamics involves 
local quenches^ where a system evolves under a Hamilto- 
nian whose eigenstates are only locally different from the 
initial state. Such evolutions (for short-ranged systems), 
as a consequence of the Lieb-Robinson bound, generate 
much less entanglement [23, 24 than those after global 
quenches, and are thus treatable with TN techniques. 

The above-mentioned Lieb-Robinson bound has 
proven essential for understanding the complexity of 
quantum states [25l [26] . It applies to systems with suf- 
ficiently short-range interactions, allowing to formulate 
a variety of general theorems, e.g., connecting excitation 
gaps and decays of correlations [27l [28]. In this article, 
we want to address the question how the Lieb-Robinson 
bound breaks down in long-range interacting systems. 

In its essence, the Lieb-Robinson bound formulates 
the principle of causality by predicting a finite speed at 
which correlations can spread. In this way, one can iden- 
tify a 'light cone' outside of which correlations are expo- 
nentially suppressed. Mathematically, under certain as- 
sumptions, the Lieb-Robinson bound can be expressed 
as a bound for the time-dependent commutator between 
two operators Oa^ Osit)^ defined at t = on two disjoint 
regions of the system A and B separated by a distance 
L |25l|26|, 
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where on the right hand side the norm is the operator 
norm, v the Lieb-Robinson velocity, and g{L) an expo- 
nentiahy decaying function. 

The Lieb-Robinson bound has important imphcations 
for thermalization: if the system locahy equiUbrates to a 
Generahzed Gibbs ensemble, also time-dependent corre- 
lation functions can be described by the same ensemble 
[29j. In some systems, the Lieb-Robinson bound can 
be understood using an intuitive pseudo-particle picture. 
This applies if the low-lying excitations can be obtained 
by populating (for translational invariant systems) dif- 
ferent pseudo-particle momentum states, with the vac- 
uum characterized by the absence of pseudo-particles. 
Then, the system responds to a local perturbation by 
emitting pseudo-particles propagating at different finite 
speeds. The fastest particles, the ones defining the causal 
cone, propagate at a speed that is often identified as the 
Lieb-Robinson velocity for that specific model. 

Much less is known about how correlations spread in 
the presence of long-range interactions, although these 
become important in many different contexts. Whenever 
one has a local model where some of the constituents 
propagate much faster than the others, one can en- 
code the presence of the fast constituents in an effec- 
tive description of the slow ones involving a non-local 
interaction. A prime example is Quantum Electrody- 
namics, describing the contact interaction between pho- 
tons propagating at the speed of light and charges. In 
the non-relativistic limit, where the charges move much 
slower than the light, the presence of photons can be 
encoded in a long-range Coulomb potential between the 
charges. Theories with long-range interactions can have 
over-extensive energies [30l[3T] and are thus strongly non- 
local. In such circumstances, one would expect that con- 
cepts like causality and the locality of quasi-particle ex- 
citations should be reconsidered. 

The purpose of this manuscript is to address this is- 
sue on a theoretical level, using complementary analyt- 
ical and numerical calculations. We find three different 
dynamical regimes, with a break-down of Lieb-Robinson 
bounds at strong long-range interactions. We are able 
to explain these regimes through the above-mentioned 
pseudo-particle picture. Finally, we discuss experimental 
regimes in trapped-ion setups where our findings can be 
observed. 

In this article, we study the out-of-equilibrium dynam- 
ics generated by a long-range interaction in the simplest 
possible scenario that can be implemented in trapped- 
ions experiments [32^, namely the long-range transverse 
Ising chain (LRTI) 
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Here, a denote the usual spin- 1/2 Pauli matrices, and we 
set the fundamental energy unit [33] and the lattice spac- 
ing to unity. We consider a finite chain of L sites with 
open boundary conditions. The parameter a, responsi- 
ble for the long range of interactions, is varied within 




FIG. 1. (Color online) (Non-) light cones, (a-c) Block 
entanglement entropy ASi = Si{t) - Si{0) from TDVP {0 = 
7r/5, L = 100). (d-f) Polarization drm = (6'f) + 1/2 from 
LSWT {0 = 7r/20). For a > 2, the excitation spreads light- 
cone like, as in the short-range model. For 2 > a > 1, there 
is no well-defined wave front, but the excitation needs a finite 
time to bridge large distances. For a < 1, the excitation 
spreads immediately over the entire system. Black dashed 
lines in (d-f) denote the maximal spin-wave group velocity. 
In (f), they practically coincide with the abscissa. 



the broad limits 3 > a > that can be realized in the 
ion setups. This allows to tune from effectively short- 
range to strong long-range physics. The parameter 6 is 
varied in the range < 6 < |, corresponding to anti- 
ferromagnetic interactions. For any a > 0, the system 
has two gapped phases, a z-polarized phase for small 0, 
and a Neel-ordered phase for values of 6 c^ 7r/2. The 
two phases are separated by a line of second-order phase 
transitions, whose universality class depends on a, as dis- 
cussed in [34|. 

Although the LRTI model does not obey the bound ^ , 
which only holds for exponentially decaying Hamiltoni- 
ans, one can still find a generalized Lieb-Robinson bound 
[25] if the power-law interactions are 'reproducing.' This 
condition, equivalent to a sufficiently fast decay, is ful- 
filled for a > 1 (see Appendix), and imposes an algebraic 
decay of correlations governed by a. 

Here, we study the effects of a on the out-of- 
equilibrium dynamics of the LRTI after a local quench. 
For this, we use as initial state the ground state IV^gs) 
of Hamiltonian ([2| at specific values of and a, and at 
time t = perturb it locally; typically IV^o) = ^l/2 IV^gs)- 
To observe the response of IV^gs) to this local perturba- 
tion, we evolve |?/^o) in time with the same Hamiltonian 
([2]). In our analysis, we employ two complementary ap- 
proaches, the quasi-exact Time Dependent Variational 
principle on MPS (TDVP) and a linear spin-wave theory 
(LSWT), see Appendix. The latter involves a higher de- 
gree of approximation, and is only valid for states with a 
sufficient degree of classical order. It has the advantage 
that it can access, with lower computational cost, larger 
times than what is possible with the TDVP. In those 
regimes where LSWT can be applied, we have checked 



that the two methods provide compatible results, show- 
ing that the time evolution they describe is essentially 
semi-classical. This agreement is reasonable, given that 
IV^o) contains a single excitation. Since the excitation 
density spreads during the evolution, the assumption of 
non-interacting quasi-particles that underlies the LSWT 
improves over time. Therefore, LSWT predictions should 
be quantitatively correct over all time scales. 

Numerical results — We now study how an excitation, 
initially localized at the chain center, spreads. For TDVP 
we choose = 7r/5, which is not accessible with LSWT 
because here magnetic long-range order is strongly re- 
duced due to a nearby quantum phase transition. In 
Fig. [l^-c, we study the spread of quantum correlations 
through the block entropy Si = —"^^pY^ogpY^ where 
pf is the n-th eigenvalue of the reduced density matrix 
pi involving the spins 1, . . . /. 

As known from [34 , in the ground state of the z- 
polarized phase of the LRTI, the long-range interactions 
cause Sl/2 cx logL for a < 2. Therefore, to isolate the 
growth of the entropy generated during the time evo- 
lution, we analyze the excess of block entropy with re- 
spect to the initial state, ASi = Si{t) — 5/(0). We have 
performed the TDVP simulations with an algorithm dis- 
cussed in the Appendix that generalizes the ones avail- 
able in the literature p4H4Q] . We consider chain sizes up 
to L = 150, and we have checked that the accuracy of 
MPS with matrix sizes x ^ 200 is sufficient. 

For LSWT, we exemplify the resulting dynamics for 
= 7r/20 (Fig. [l]i-f), where the ground state is almost 
perfectly polarized, (S^) ~ —1/2 [34]. In this case, a 
useful measure for the spread of the perturbation is the 
excess magnetization drrii = (Sf) + 1/2. Notably, within 
LSWT, this gives a direct measure for the single-site en- 
tanglement entropy, S^ = {Snii -\- 1) \og{Smi + 1) — 
Srrii log 6mi [HI [42] . Figure IT] evidences the similar be- 
havior for the two methods and the two regimes. 

For generic ^, we identify three dynamical regimes as 
a function of a. (i) For a > 2 [realized in Nature, e.g., 
for van-der-Waals {a = 6) or dipole-dipole {a = 3) in- 
teractions] , the system behaves as if short-range interact- 
ing, with the spread of the local perturbation essentially 
bounded by a light-cone. The excitation maximum de- 
fines a clear wave front, and its linear propagation gives 
a constant Lieb-Robinson velocity as in the short-range 
model, coinciding with the maximal spin-wave group ve- 
locity. Outside the light cone, correlations decay alge- 
braically with a power determined by a, thus obeying 
the generalized Lieb-Robinson bounds, (ii) In the range 
2 > a > 1, although at short times there appears an 
effect resembling a light cone, it does not really bound 
the propagation of the perturbation, since correlations 
consistently leak out of it. Further, the propagation of 
the perturbation shows complex interference effects due 
to longer-range spin fiips, and at sufficiently large times 
one cannot identify an unambiguous wave front. Still 
the excitation needs a finite time to bridge larger dis- 
tances, (iii) For a < 1 (typically, a = 1 corresponds 
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FIG. 2. (Color online) (a) Spin-Avave dispersion rela- 
tions. Inset: For a > 2, cjfc is a deformed cosine, similar to 
the short-range case. For a < 2, cjfc develops a cusp at /c = tt, 
which becomes sharper for smaller a. Main panel: For 
a < 1, the number of modes with |t'^| = \duJk/dk\ > {L/2)/to 
diverges with L for any to > 0. In the picture, to = 50; cir- 
cles: L = 20 (6 modes); triangles: L = 40 (8 modes), (b) 
Maximal group velocity for different L (at ^ = 7r/20). 
Top: For a > 2, '^max is essentially constant as a function of 
a. At a < 2, t'max increases sharply. Bottom: For a > 1, 
Vmax/L ^ for L ^ oo. The time tb = L/{2vmax) at which 
excitations reach the system boundary diverges. For a < 1, 
Vmax/L increases with system size. Information about the 
local quench reaches the entire system instantaneously. 



to Coulomb- or gravitation-like potentials), the general- 
ized Lieb-Robinson bounds valid for a > 1 can no longer 
be defined. As a consequence, the system becomes truly 
long ranged, and correlations spread practically instan- 
taneously over the entire chain. 

These results complement the one available about ther- 
malization in disordered systems with random interac- 
tions that are modulated by a long-range power law. In 
that case, the time average of local observables tends to 
a value predicted by a Generalized Gibbs ensemble only 
if a< 1 [43. 

Our findings differ from previous results for the spe- 
cific cases of Hamiltonians made of mutually commuting 
terms, such as, e.g., Eq. Q with = 7r/2. In such set- 
tings, at a = 1 nothing special seems to happen, while the 
value a = 0.5 seems to separate two dynamical regimes 
^44j, in one of which one finds, e.g., prethermalization 
plateaus [45]. Also, during time evolution with commut- 
ing Hamiltonians, the block entropy of subsystems can 
increase unchecked with block size for a < 0.5, whereas 
for a > 1 it is strictly upper bounded [46]. 

Pseudo-particle dispersion relation — The qualitatively 
different behavior in the regimes (i-iii) can be under- 
stood in a simple quasi-particle picture: During the lo- 
cal quench, all spin-wave /c-modes become approximately 
equally populated (occupation ^ ^/L). If the pseudo- 
particles do not interact (a good approximation for low 
pseudo-particle density), they will subsequently propa- 
gate with the group velocity corresponding to their k 
value, Vg = ^^, which depends only on the dispersion 
relation ujk (see Fig. [2k; see Appendix for an analytical 
formula from LSWT). 

In the range 2 < a < oo, the maximal group veloc- 
ity Vmax is achieved around /c = tt, and does barely de- 



pend on system size or a (Fig. [2b, top). At a < 2, 



however, uok acquires a cusp at /c = tt, where 
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verges hke |/c — 7r| ~ . Consequently, v^ax is attained at 
/c = TT =b 27r/I/ [47], and diverges as v^ax oc (27r/L)"~^. 
Stih, the time scale in which pseudo-particles can reach 
the boundary, tb = L/{2vmax)^ scales as L"~^, which di- 
verges for 1 < Of < 2; the time to reach the boundary 
increases with system size, even for the fastest mode. 
The long-range effects become more dramatic at a < 1, 

with a stronger divergence Vg oc^ |/c — 7r|^"~ ^ . Here, 
for the fastest mode, t^ decreases with system size (in 
fact, even for a diverging number of modes, see Fig. [2] 
and the Appendix). In Fig. [2J3, the transition between 
the three regimes can be clearly identified. 

The spin- wave dispersion also explains the diffusive ef- 
fect encountered at small a (see Fig.[lfe-f). With decreas- 
ing a, the dispersion becomes fiatter around the sides of 
the Brillouin zone. Therefore, there are many slow quasi- 
particles that remain in the central region for a long time, 
giving rise to an apparent diffusive core of high density. 

Scaling of entanglement entropy — To numerically con- 
firm the validity of the pseudo-particle picture, we ana- 
lyze within the TDVP the scaling of the EE of half of 
the chain SL/2{t) during the time evolution. Interest- 
ingly, for all values of a considered, the excess entropy 
ASL/2{t) initially increases as a power of t and then sat- 
urates to a value very close to ASL/2{t) = log 2, inde- 
pendent of system size (Fig. [sk). The initial growth is 
faster for smaller a, in agreement with the presence of 
faster pseudo-particles. Before entering the saturation 
regime, however, systems with smaller a start to evolve 
slower, in agreement with the appearance of a diffusive 
evolution. The fact that the excess of EE of a block 
saturates to a value independent of its size is in remark- 
able contrast to the ground-state properties. This effect 
finds a natural explanation in the semi-classical picture 
of pseudo-particles: the states that dominate the time 
evolution are states with only one pseudo-particle; the 
log 2 is then immediately understood as coming from the 
two orthogonal possibilities of the pseudo-particle being 
either in the left or in the right half-chain. 

A further confirmation comes from the half- 
chain entanglement-spectrum evolution, hn{L/2,t) = 
logp2/2(^)5 where p2/2 ^^ ^^^ n-th eigenvalue of the re- 
duced density matrix of half of the chain. As seen in 
Fig. [sJd, the spectrum is dominated by only few eigen- 
values, with two of order one as expected from the log 2 
asymptote, and a huge number of eigenvalues below 10~^. 
These eigenvalues grow steadily, but we expect that they 
do not affect equilibrium properties, since they are asso- 
ciated to higher energies and thus, at long times, their ef- 
fect should average out. These findings are in agreement 
with similar observations in short-range systems [24l [48] , 
where semi-classical models provided a good description 
of these kinds of out-of-equilibrium dynamics. 

Experimental implementation — Finally, let us remark 
that our system does not break non-locality of the phys- 
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FIG. 3. (Color online) (a) Growth of the entanglement 

entropy. The excess ASL/2{t) grows initially as a power law 
with t for all considered a. It then saturates to log 2 indepen- 
dent of system size, as expected from the pseudo-particle pic- 
ture. This is shown in the insets where we compare the satu- 
ration value for chains of different length and for completeness 
show that there is no residual dependence of the saturation 
value on the MPS matrix dimension x- (t>) Evolution of 
the entanglement spectrum. The entanglement spectrum 
is dominated by only two eigenvalues, which in the pseudo- 
particle picture correspond to the pseudo-particle being in 
the left or the right part of the chain. The other eigenvalues 
are significantly smaller, as a confirmation of the goodness of 
semi-classical descriptions of the evolution. 



ical world, as one should hope. E.g., in the trapped- 
ion implementation, Hamiltonian ^ describes an effec- 
tive dynamics for electronic states of the ions, which are 
coupled by collective phonon modes by employing laser 
fields [32l [49l ISO] . The phonon dynamics can be neglected 
on time scales much larger than those associated to the 
detuning between laser driving and phonon frequencies. 
These time scales are typically O{10fis). Moreover, the 
derivation of Eq. ^ employs a rotating-wave approxi- 
mation in the phonon frequencies, corresponding to ne- 
glecting terms that average to zero on time scales (9(l/is). 
When the group velocity reaches these time scales, the 
effective Hamiltonian S breaks down, just as how the 
Coulomb potential is no longer valid when charged par- 
ticles move close to the speed of light. On the other 
hand, the time scale of the spin interactions is typically 
h/J = O(lms). Therefore, although the group veloci- 
ties of the spin system cannot truly diverge, they can be 
several times larger than the scale set by J, thus still 
providing a drastic effect. Therefore, in typical practical 
implementations [9-14 , it should be possible to explore 
the three widely different regimes discussed above. 

In typical experiments, the system size is finite, but the 
three regimes are distinguishable already in small chains. 
In the short-range regime (i), the ratio of wave- front max- 
imum and the subsequent minimum increases with sys- 
tem size, thus defining an increasingly sharp wave front. 
Contrarily, in the weak long-range regime (ii), this ratio 
decreases with system size, until the wave- front maxi- 
mum disappears, making it impossible to unambiguously 
define a wave front. In the strong long-range regime (iii), 
the excitations start to reach the boundary of the system 



immediately, independent of system size. 

Conclusions — We have numerically identified three 
qualitatively different regimes of the dynamics in the 
LRTI model, and we have reached an intuitive under- 
standing through analytical studies of the spin- wave dis- 
persion relation. The validity of such studies is backed 
by the agreement of LSWT results with results extracted 
with the quasi-exact TDVP in the short- time regime. 

We have found qualitative changes of the spin-wave 
dispersion relation at a = 2 and a = 1, indicating 
the transition from short-range to weakly long-range to 
strong long-range physics. In the last case, a strong diver- 
gence of the quasi-particle velocities leads to a practically 
instantaneous spread of excitations through the system. 
It will be extremely interesting to study how this effect 
influences thermalization and how our findings carry over 
to more complicated models or other dimensionalities. 

Finally, we have discussed the regime where strong 
long-range physics should appear in typical trapped-ion 
experiments, and outlined a way to observe it in small 
chains. We hope our findings to inspire experiments 
along these lines. 
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Marie Curie project FP7-PEOPLE-2010-IIF ENGAGES 
273524, TOQATA (FIS2008-00784), Spanish MICINN 
(FIS2008-00784), Catalunya-Caixa, the Austrian Science 
Fund (SFB F40 FOQUS), and the DARPA OLE pro- 
gram. 



Appendix A: Time-reversal-symmetric scheme for 
the TDVP applied to long-range Hamiltonians 

In a previous work, it has been shown how to find an 
MPS approximation to the ground state of a long-range 
Hamiltonian using an extension of the time-dependent 
variational principle (TDVP) [34 . The TDVP is an algo- 
rithm that uses the geometric notion of the MPS tangent 
plane [38l [40] . It allows to find the ground-state descrip- 
tion as a MPS by solving a differential equation for the 
tensors defining the MPS. The same generalization of the 
TDVP presented in [34j can be used to perform real-time 
evolutions [38l |40] for systems with long-range interac- 
tion, which we have exploited to compute the quench dy- 
namics described in the main text. The main difference 
when performing the dynamics - instead of imaginary 
time as is necessary to obtain the ground state - in real 
time, is that special care has to be taken to ensure that 
the algorithm does not violate time-reversal symmetry. 
This immediately ensures that the algorithm conserves 
both the norm and the energy of the initial state, as it 
should. Care is advised, as real-time evolution does not 
enjoy the same self-correction as imaginary-time evolu- 
tions where small errors in one step can be corrected in 
the next step. Here, instead, any small error is propa- 
gated along the simulation, and one needs further care 



to minimize those. 

Here, we describe the steps we follow in order to ensure 
the time-reversal invariance of the integrators scheme 
that we apply to the above-mentioned differential equa- 
tions. This is a simple modification of the technique pro- 
posed in [38l [51] , suitable for evolutions dictated by long- 
range Hamiltonians in finite chains. 

First, we briefly recall the general strategy that is com- 
mon to both real-time and imaginary-time dynamics. 

(1) Encode the starting state of the time evolution 
as a MPS described by the set of tensors L {An}, 
An^n = 1---L |'0o{^}) and (2) encode the long-range 
Hamiltonian as a MPO described by a set of tensors L 
{O}, H{{0}). The evolved state is obtained by (3) solv- 
ing the Schrodinger equation for a short time interval dt 
with initial condition given by IV^o), 



idtm{A{t)})) = H{{o})m{A{t)})), 



(Al) 



where we set h = 1. (4) To solve the above equation in- 
side the manifold of MPS with fixed tensor dimensions, 
one needs to introduce tangent vectors. These are gener- 
ically defined through two sets of tensors {A} and { 5 }, 
and are expressed as the linear combination of MPS de- 
fined by A everywhere but one Bn at a specific place, in 
formula 

\T{{A} ,{B})) = J2m{AU,B„)) , (A2) 

n 

where we have used the notation { ^ }^ to define the set 
of A tensors from which we have removed tensor A^ . 



(5) While the left hand side of the above equation ( Al ) 
defines a tangent space to the manifold of the MPS with 
fixed bond dimension, the right hand side is not contained 
in that space and should be explicitly projected onto it. 
In formula, we would like to find the tangent vector \T) 
that minimizes the distance from H \ip{{A{t)})), 



\T{A},{B* 



}) , : min I 



\T)-Hm{A{t)}))\f. (A3) 



In practice, in the canonical form, the computation is 
simplified by requiring that the tangent vectors are or- 
thogonal to the original vector. To ensure the orthogo- 
nality, the Bn tensors in the tangent vectors are defined 
as the contraction of auxiliary tensors, (for normaliza- 
tion convenience) the inverse square root of the reduced 
density matrix, times a matrix of free coefficients of di- 
mension called X^, and a fixed projector Vn on the or- 
thogonal space to the one on which the starting vector is 
defined. 



(6) At this point, one can discretize Eq. (Al) and in- 
tegrate it iteratively through 



An{t^dt)=An{t)^idt{Bny, 



(A4) 



for all sites n = 1 . . . L. 

We now turn to real-time dynamics, and we focus 
specifically on designing a time-reversal-invariant inte- 
grator scheme. This requires improving the first-order 



integrator ( |A4[ ) to at least the so called middle-point in- 
tegrator. This involves finding an intermediate step for 
each n, An{t-\- dt/2) such that both An{t-\-dt) = An{t-\- 
dt/2)Mdt/2{BnY and An (t) = An{t^dt/2)-idt/2{BnY , 
where { 5 } is the set of tensors defining the projection 
onto the tangent space of the action of the Hamiltonian 
on the state \ilj{{A(t + dt/2) })). The two above condi- 
tions can be taken as a definition of An{t^dt/2) that we 
then use to complete the evolution step by just integrat- 
ing the state for another dt/2^ 



An{t + dt) = An(t + dt/2) + idt/2{Bny, 



(A5) 



so that we ensure that the evolution is invariant under 
time reversal. The important part becomes finding the 
intermediate An{t -\- dt/2). As suggested in [38 , one can 
devise an iterative procedure to determine An{t + dt/2). 
Here, we describe an alternative procedure to the one pre- 
sented in [38 that is well suited for finite-chain Hamil- 
tonians encoded in MPO as the ones discussed in this 
paper. The procedure consists in proceeding locally in 
the evolution (site by site) requiring that each local step 
is time-reversal invariant. For this reason, chosen a posi- 
tion n in the chain, one proceeds by 

(1) Obtaining a trial A^{t + dt/2) by solving Eq. ( [AT ) 
for a time step dt/2, with the initial state locally de- 
scribed by An{t). (2) Obtaining a trial B by finding 
the best tangent vector that approximates the r.h.s. of 
Eq. ([Al|. (3) Evolving back A^(t + dt/2) to A^t) 



(that initially will differ from An{t)) by solving Eq. (Al) 
for a time step —dt/2, with an initial state locally de- 
scribed by An{t + dt/2). (4) Compute the AA^{t) = 
Ai{t) — A^{t) and project it onto the tangent space de- 
fined at A^{t + dt/2). (5) Compute the error as E^ = 

V\\ IV^({^}i'^^?(^))) II- (^) I^ ^his way, we can obtain 
the improved estimate of the middle point Ai{t + dt/2) 
as A}{t + dt/2) = A^(t + dt/2) + P^o AA^(t), where P^o 
is the projection on the tangent plane at the old estimate 
A°^{t + dt/2). 

We then repeat the procedure starting again from step 
(2) and iterate as often as necessary in order to bring the 
error in the inversion En below the required precision 



(typically around 10~^^). The procedure is repeated for 
all sites, and at the end of a sweep from 1 to L, one 
completes an elementary evolution step of dt. For more 
details about the other aspects of the algorithm and pos- 
sible improvement using higher-order Ruge-Kutta inte- 
gration schemes, we refer the reader to the literature on 
the subject [34| [38l [40] . 



Appendix B: Linear spin- wave theory for long-range 
models 



To gain some analytical understanding of the dynam- 
ics of the long-range system described by Hamiltonian 
([2|, we employ a linear spin- wave theory (LSWT). This 
theory is well known to yield good qualitative results 
in phases with strong magnetic order [52], such as the 
strongly z-polarized phase occurring for small [3T. In 
our numerical analysis, we will therefore focus on that 
case, although we will keep our derivations general. An 
advantage of LSWT is that the long-range interactions 
are implemented into the formalism without additional 
complications, as we will sketch now. 



1. Determining the ground state 

As a first step to finding the ground state of spin waves, 
it is convenient to rotate the spins into a local, twisted 



coordinate system 



x',y', 



so that the new 



z axis IS 



aligned with the quantization axis. In the antiferromag- 
netic case of > 0, a convenient form is to rotate the 
spin 1/2 operators into S^ =lZiSi, where 



Ui 



( — 1)* COS7 


(-1)^+1 sin 7 








- sm7 

cos 7 



(Al) 



Since only S^ and S^ operators occur in Hamiltonian ([2| , 
it is sufficient to restrict the rotation to the xy plane. We 
keep the angle 7 free at this stage and will find it later 
through the minimum of the energy. 



In terms of the rotated spin operators, the system Hamiltonian reads 



i:^ = 4sin6>^ 



{ij) 



\^-J^ 



(-1)^+^' cos^-fSfSf + (-l)^+^'+i sin7cos7(5f ^f + S^' Sf) + (-1)^+^' sin" -fSf S', 



2 cos 6> ^ (sin -iSf + cos -^Sf) . 



(A2) 



For a state that is strongly polarized along the z' axis, 
one can approximate spin-S* operators (here S = 1/2) by 
bosonic operators via the Holstein-Primakoff transfor- 



mation 52 , S'f ^ S'- 



-a]ai 



5+ 



-^ 



1 



2S ' 



and 



S- ^V2S\l-^ 



2S 



We now insert these into Hamil- 



tonian (A2), neglect contributions beyond linear order in 



^^, and use that terms that are linear in the boson op- 
erators vanish in the minimum of the free energy. More- 



over, we apply a Fourier transform a J = -j^^k^ 



:^ikri^\ 



fc' 



leading finally to 



H = 2_.\%^k 2{2Ssm0cos^ Ilk + cos 6> cos 7 — 46'sin6>sin^ 77Q ^) + {al.a_j^ + a/ca_/c)25'sin^cos^ 77^"^ 

k 

+2^ sin cos^ 7 Y^ 7^"^ + L{2Sf sin sin^ 7 ^^^""^ - L2S cos 6> cos 7 , 



(A3) 



where we defined 



t'=T.^-^^oskS. 



(A4) 



6>0 



This last abbreviation encorporates the entire long-range 
nature of the system, preserving the extreme simplicity 



I 

and elegance of LSWT. 

Hamiltonian (A3) can now be diagonalized as usual 
by a Bogolioubov transformation, ak = cosh /3/e a^ + 

sinh/3/c aL/g, cl_]^ = sinh/3/c a/c + cosh/3/c aL/.. Demand- 
ing that the a^ obey bosonic commutation relations, and 
that only terms proportional to a^at yield a contribu- 
tion to H^ one obtains the Bogolioubov angles cosh 2/3/e = 
Bk/uJk, sinh2/3/e = -2Ak/uJk, and 



H = Y,^k UW + ^) + IZ (^fc - \Bk - 25cos^cos7 + (25)'sin^sin2 77^"M , 



(A5) 



where 



B]^ = 4S sin cos^ 7 7^"^ + 2 cos cos 7 — 85 sin sin^ 7 7q ^ , 

Ak = 2S'sini9cos^7 7^''\ 



(A6a) 
(A6b) 



and with dispersion relation 



^k 



Bl-AAl 



(A7) 



The ground state of Hamiltonian (A5), I'^gs)^ is found as 
the vacuum of Bogolioubov particles, ak\^) =0 \/k. We 
can now determine the free-energy minimum 7 by min- 
imizing ('0Gs| ^ IV^Gs)- Alternatively, one can demand 
that terms that are linear in the spin- wave operators van- 
ish at the energy minimum, which gives a condition on 7 
as a function of 0. For the case of ^ = 7r/20 that we use 
in the main text, we find 7 = independent of a, and 
the spins are strongly polarized in negative z direction. 



2. Spin-Avave group velocity 



The dispersion relation (A7) determines the group ve- 



locity, which is of the form 

dujk 






^2% 



(A8) 



C3 



with Ci^2,3 constants. Divergences of Vg (as given in the 
main text, see Fig. ^ can hence be found easily by ana- 



lyzing 7^"^ and 




dk 


.5>0 



sin k5 . 



(A9) 



As illustrated in the main text for the example of 
9 = 7r/20 (Fig. ^), there are two transitions that can 
be found generically through an analysis of ^max = 
max/e'u^(A:) as a function of a. For a > 2, 'Umax is al- 
most constant as a function of a, whereas below it, it 
rises rather steeply with decreasing a. This indicates a 
transition in the dynamical behavior of the system. Ad- 
ditionally, for a ^ 2 the /c-value where Vmax is achieved 
lies around k = 7r/2 and changes slowly with decreasing 
Of, whereas at a = 2, it transitions to /c = tt (Fig. |Al[ left 
panel). 

In the range 1 < a < 2, although Vmax achieves large 
values, Vmax/L scales to zero with increasing system size, 
indicating the locality of information in these systems. 
This changes drastically at a < 1, where Vma^/L in- 
creases with system size (see discussion in the main text). 
In that regime, the fastest mode reaches the boundaries 
earlier for larger systems; the information is distributed 
essentially instantaneously over the entire chain. Actu- 
ally, the number of modes diverges for which the time to 
reach the boundary decreases with system size. Conse- 



k(v™'')/jt 




3. Linear spin-w^ave theory and dynamics 



To evaluate the time evolution under H^ we make use 
of the fact that all involved states remain Gaussian at all 



times. Since Hamiltonian (A3) is quadratic in the boson 



operators, its ground state is completely determined by 
the correlators 



FIG. Al. (Golor online) Properties of the maximal 
group velocity. Left: The k value at which t'max is found 
makes a transition to k — n dX a — 2. Data is for L — 500. 
Right: Where the constant ci vanishes, the system is dis- 
persionless. This happens at the points = 0,7r/2,7r, where 
there are no non-commuting interactions in the Hamiltonian. 
Then, the system dynamics is localized. 



quently, Vto > 0, there exists a chain length Lq such that 
ML > Lq the number of /c-modes with Vg{k) > {L/2)/to 
is larger than any given no = cL^^~"^/^^~"\ with c a 
constant. For a < 1, no diverges with L. In other words, 
we find a diverging number of quasi-particle modes that 
reach the boundary before the (arbitrarily small) time 
to. The physical reason is simple: the mode spac- 
ing decreases faster with L than the /c-interval where 
Vg{k) > {L/2)/to (see symbols in the main panel of 
Fig#). 

This behavior is generically true, except when in 
Eq. (A8) we have Ci = SS sin cos ^ 7( cos i 
AS sin sin^ 77o 



Al 



cos 7 — 
right panel. 



) = 0. As seen in Fig. 
this happens only at ^ = 0, 7r/2, tt, . . . . At these pa- 
rameter values, one has either only the magnetic field or 
only the spin-spin interactions, i.e., there are no non- 
commuting terms in the Hamiltonian. In that case, inde- 

pendent of divergences of ^^ , the quasi-particle disper- 
sion relation ( |A7[ ) becomes dispersionless, and the system 
dynamics is localized. 



Z ZIj ^ UJk 

(AlOa) 



Gij = (V^Gsl «1«] IV^Gs) = lyyY^ 



-2A, , 



k{rj—ri) 



2L ^ ujk 



(AlOb) 



At time t = 0, we quench the system with 
a spin fiip at site tti, corresponding in LSWT 
to the operator S!^ = ^{am + ciln)^ ^^^ ^^^ 
state becomes |?/^o) = S!^ IV^gs) /vW, where 

JV= [3? ((V^GS I amttm I V^GS)) + (V^GS I aJn^m I V^GS) + ^] /2 

is the normalization [53 . Since the initial state is 
Gaussian, expectation values after the quench such 
as {alaj){t = 0) = (?Ags| 5'^'4%^m IV^Gs) /A/" can be 
decomposed into a combination of expectation values of 
two-point correlations before the quench, using Wick's 
theorem [54]. For example, (V^gsI <^Jn^I%^m IV^Gs) = 

(V^GS I «]n<^i IV^GS) (V^GSI ^ttJn iV^Gs) + 

(V^GSI ttJn^ IV^GS) (V^GSI ^^L iV^Gs) + 

(V^Gsl «]n«L iV^Gs) (^Gsl CLJaj H^Qs) - Since the corre- 
lations after the quench remain Gaussian, and since a 
Gaussian state remains Gaussian under the application 
of a quadratic Hamiltonian, it is sufficient to consider 
the time evolution of the two-point correlators. 



To compute the time evolution, we use Heisenberg's equation of motion, which for an arbitrary operator A reads 
in the small-time limit, A{t + At) = A{t) + iAt [H, A]. We find 



Fab{t + At) =Fab + iAt2^sin^cos^7[^ . J^ {G^ + F,,) - ^ 



(-1) 



i-b 



\i-b\' 

i-iy 



(G. 



a V ci* n^ -^ o2 



F„, 



(Alia) 



Gab{t + At) =Gab + iAt{2SsmeC0S^ 7 [^ . _' ^ {Gib + Fbi) + ^ . J « {Gia + Fai) 

+ 4 (cos d cos 7 — 45 sin d sin^ 7 7q ^ ) Gab \ 



(Allb) 



where the right hand side is to be evaluated at time t. In our numerical evaluation, we chose At = 0.002, and checked 
that a further decrease does not improve the results. In the main text, we plot as a function of time the deviation of 
the magnetization from —1/2, ^vni = (5'f ) + 1/2, which is nothing else than Fait) — 1/2. 

I 



Appendix C: Power-law interactions are reproducing 
if and only if a > 1 



interactions that decay at least exponentially with dis- 
tance 1—3 between lattice sites i and j. For interactions 



Typical Lieb-Robinson bounds with the associated ve- 
locity are defined for short-range interacting systems, i.e.. 



K{i — j) that decay slower than exponential, one can 
still define Lieb-Robinson bounds on the commutators 
of operators — similar to the bound (IT]), but without a 
constant Lieb-Robinson velocity — as long as K{i — j) 
is reproducing [25,, i.e., if it fulfills the condition 



E 



,my^ij 



K{i — m)K{m — j) < XK{i — j) \/i^j 

(Al) 

for some constant A. For simplicity, we discuss here open 
boundary condition for a chain of length L -\- 1 with L 
even. For a power law K(i — j) = 1/ |i — j|", this con- 
dition is fulfilled if the decay is faster than a > 1 and 
violated for a < 1, as we will show now. 



1. Power- law interactions are reproducing if a > 1 



It is convenient to rewrite the condition ( Al ) using the 
definition 



P{iJ) 






■Jl 



■ to| I'm — j\ 



(A2) 



SO that it becomes P{i^j) < A. To show that 1/ \i — j\^ 
is reproducing for a > 1, we need to demonstrate that 
P{i^j) converges with L for any z, j. Consider i and j 
placed symmetrically at positions ±(^/2 for S even. Then, 



S S 
2'2 



E 

\m\> 



(5" 



, (m+f)"(m-|)" 

2 

(5" 



iJ^,(f + -)"(! 



(|+m)"(| -m)" 



(A3) 



Let us treat the two sums separately. The first sum 
is upper bounded by E|mi>^/2 (m.S/9.^2c. ^ the last term 



^\m\>6/2 (m-5/2)2« 



of which reads M 



4«5« 



For constant (5, this goes 



(L-5)2«- 

to zero as ex L~^". From this, it would seem that this 
sum converges for a > 1/2. However, one can consider a 
more demanding scenario, which shows that convergence 
is reached only for a > 1, namely, if one lets 5 increase 
with system size, 6 = d{L) = cL^ . Here, the condition 
S < L demands /3 < 1. Then, 



M: 



^a^aj^alS 



< 



4"c^ 



L2«(l - cL/5-l)2« - la^l _ ^)2a ' 



(A4) 



where we used f3 < 1 to bound L"^ < L" and cL^~^ < e 
(with e arbitrarily small for /3 < 1, provided L is suffi- 
ciently large, and e = c < 1 for /3 = 1). Therefore, M 
decays at least as fast as L~", meaning that the sum over 
it is assured to converge for a > 1. 



For constant ^, the second sum in Eq. (A3) is constant. 



The only way it can increase is by increasing 5 in some 
way with L. To study if this can make it diverge, consider 
the difference of when it is evaluated at S and S -\- 2, 



E 



{s+2r 



, , , (| + l + m)"(| + l-m)^ 

|m|<|+l ^2 ' ^2 I 



(5" 



(^ 

0<TO<| ^2 



\m\< _ 

=2 E 

0<m< 
(1 + 1)2" 

+ E 



(|+m)«(| 
1 






((^ + 2)^ 
(f +m + 2)^ 



5^ 



(f +m)« 



(A5) 



((^ + 2)« 



^(f + l + m)«(f + l-m)« ' 



where we regrouped the terms of the two sums into con- 
tributions from the outermost summands, the one at the 
origin, and the additional summands that are inserted 
beside the origin upon increasing 5. One can show for 
the first sum of this expression that all terms are neg- 
ative, as is the case for the contribution at the origin. 
Now, we only have to show that the last few terms de- 
cay sufficiently fast. In fact, they decay as (5~", so that 
even increasing 5 proportional to L leads to a convergent 
sum as long as a > 1. Therefore, in this case, also the 
second sum in Eq. ( |A3| ) converges. The argumentation 
here can be carried over to positions deviating from the 
symmetric case i, j = ±(5/2. We have thus demonstrated 
that K{i — j) (X l/|z — j|" is reproducing for a > 1. 



2. Power-law interactions are non- reproducing if 

a < 1 

To show that K{i — j) is non-reproducing for a < 1, it 
is sufficient to demonstrate the divergence of P{i^j) for 
a specific case, which can easily be done for the choice 
i = -L/2, j = L/2, 



L/2-1 

P{-L/2,L/2)= Yl 

m=-L/2+l 
L-1 



L^ 



(m + L/2)«(m-L/2)" 

I L-l ^ 

The last sum converges towards the Riemann-zeta func- 
tion C{a). This lower bound for K{i — j), therefore, 
diverges for a < 1, where K(i — j) is hence non- 
reproducing. This property explains the violation of the 
Lieb-Robinson bounds for small a. 
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